{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00000-eb52ae5e-a1f6-4ad2-99c6-d6eec6b07cfe"
   },
   "source": [
    "# 3D dose map analysis\n",
    "\n",
    "Use the output of the following simulation:\n",
    "- Folder: dosimetry/\n",
    "- Macros: ex2.mac, ex3.mac\n",
    "\n",
    "Helping ressources: http://insightsoftwareconsortium.github.io/SimpleITK-Notebooks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "cell_id": "00001-6d57c54d-e209-4062-a0b1-7e0992fe9b6f",
    "execution_millis": 1,
    "execution_start": 1605009963671,
    "output_cleared": false,
    "source_hash": "4786f891"
   },
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import SimpleITK as sitk\n",
    "from ipywidgets import interact\n",
    "import os\n",
    "from pathlib import Path"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to be able to run the interactive plot, you must install ipympl as indicated here: https://github.com/matplotlib/ipympl\n",
    "This is activated by the previous `%matplotlib widget` line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "cell_id": "00002-fee9b460-446b-4c59-ad8e-e624b5c380a8",
    "execution_millis": 1,
    "execution_start": 1605010049758,
    "output_cleared": false,
    "source_hash": "702ff6d0"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Current Working Directory (CWD) is: \n",
      " /Users/dsarrut/src/gate/dqprm/2021/gate-exercices/dosimetry\n"
     ]
    }
   ],
   "source": [
    "# The following command display the current working directory (where jupyter has been launched)\n",
    "cwd = os.getcwd()\n",
    "print('The Current Working Directory (CWD) is: \\n', cwd)\n",
    "folder = Path()\n",
    "if (not folder.is_dir()):\n",
    "    print('ERROR: {} is not a folder.'.format(folder))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "cell_id": "00003-218500bc-1c37-4681-8413-9feff467bf2b",
    "execution_millis": 0,
    "execution_start": 1605010082790,
    "output_cleared": false,
    "source_hash": "77b2d837"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image size:  (126, 126, 111)\n",
      "Image spacing:  (2.0, 2.0, 2.0)\n",
      "Image origin:  (0.0, 0.0, 0.0)\n"
     ]
    }
   ],
   "source": [
    "# Read a sitk image\n",
    "filename = os.path.join(folder,'./data/patient-2mm.mhd')\n",
    "img_ct = sitk.ReadImage(filename)\n",
    "print('Image size: ', img_ct.GetSize())\n",
    "print('Image spacing: ', img_ct.GetSpacing())\n",
    "print('Image origin: ', img_ct.GetOrigin())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "cell_id": "00004-2d061556-d95d-42a1-8970-04f5211dc420",
    "execution_millis": 4,
    "execution_start": 1605010084508,
    "output_cleared": false,
    "source_hash": "672e2d66"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Array size:  (111, 126, 126)  <--- be careful to the dimension order!)\n"
     ]
    }
   ],
   "source": [
    "# Convert sitk image to a numpy array\n",
    "arr_ct = sitk.GetArrayFromImage(img_ct)\n",
    "print('Array size: ', arr_ct.shape, ' <--- be careful to the dimension order!)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "cell_id": "00005-d092e942-514e-4700-95ea-1be0f120dfed",
    "execution_millis": 302,
    "execution_start": 1605010085960,
    "output_cleared": false,
    "source_hash": "1a513e72"
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1a78fff025de45f59a4fea3b01685528",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "287620e6713e46c89ec7a5a61ab42cc3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=63, description='sx', max=126), IntSlider(value=63, description='sy', ma…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,3, figsize=(8, 2))\n",
    "def show_ct(sx,sy,sz):\n",
    "    ax[0].imshow(arr_ct[sz,:,:], cmap=plt.cm.gray)\n",
    "    ax[1].imshow(arr_ct[:,sx,:], cmap=plt.cm.gray)\n",
    "    ax[2].imshow(arr_ct[:,:,sy], cmap=plt.cm.gray)\n",
    "    \n",
    "interact(show_ct, sx=(0,img_ct.GetSize()[0]), sy=(0,img_ct.GetSize()[1]), sz=(0,img_ct.GetSize()[2]));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "cell_id": "00006-aedcaa29-ceb7-4460-9e45-b3b141ee79b1",
    "execution_millis": 4,
    "execution_start": 1605010087540,
    "output_cleared": false,
    "source_hash": "7261243f"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image size =  (50, 50, 50)\n",
      "Image min and max:  0.0 2.5758645e-05\n"
     ]
    }
   ],
   "source": [
    "filename = os.path.join(folder, './output_ref/3d-pat-proton-Dose.mhd')    # <--------------------------------------TO CHANGE BY YOUR OWN OUTPUT FOLDER \n",
    "img_dose = sitk.ReadImage(filename)\n",
    "arr_dose = sitk.GetArrayFromImage(img_dose)\n",
    "print('Image size = ', arr_dose.shape)\n",
    "print('Image min and max: ',  np.amin(arr_dose), np.amax(arr_dose))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "cell_id": "00007-2cbb61a8-ddc8-4365-b8d9-68d5f4d0656e",
    "execution_millis": 2,
    "execution_start": 1605010088048,
    "output_cleared": false,
    "source_hash": "98f4f536"
   },
   "outputs": [],
   "source": [
    "filter = sitk.RescaleIntensityImageFilter()\n",
    "filter.SetOutputMaximum(1.0)\n",
    "filter.SetOutputMinimum(0.0)\n",
    "img_dose = filter.Execute(img_dose)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "cell_id": "00008-c500cafb-45e4-42ea-bd48-71b432ecb0b8",
    "execution_millis": 2,
    "execution_start": 1605010089021,
    "output_cleared": false,
    "source_hash": "29c4ff45"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image size =  (50, 50, 50)\n",
      "Image min and max:  0.0 1.0\n"
     ]
    }
   ],
   "source": [
    "arr_dose = sitk.GetArrayFromImage(img_dose)\n",
    "print('Image size = ', arr_dose.shape)\n",
    "print('Image min and max: ',  np.amin(arr_dose), np.amax(arr_dose))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "cell_id": "00009-c0d6b0ff-c008-4dcb-8ad6-b9c0cd7f68c0",
    "execution_millis": 167,
    "execution_start": 1605009991135,
    "output_cleared": false,
    "source_hash": "6e7db63d"
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c8e30318b0bf44eb9f1e9d7857c84869",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1e0132417e6c48c9babc659cd8d98abe",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=24, description='nslice', max=49), Output()), _dom_classes=('widget-inte…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig1, ax1 = plt.subplots()\n",
    "def show_dose(nslice):\n",
    "    im = ax1.imshow(arr_dose[:,nslice,:], cmap=plt.cm.gray)\n",
    "interact(show_dose, nslice=(0,len(arr_dose)-1));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "cell_id": "00010-31e68897-2f11-40f2-8942-eaccc9e003dd",
    "execution_millis": 11,
    "execution_start": 1605009994568,
    "output_cleared": false,
    "source_hash": "7aba1c68"
   },
   "outputs": [],
   "source": [
    "img_resampled_dose = sitk.Resample(img_dose, img_ct, sitk.Transform(), sitk.sitkLinear, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "cell_id": "00011-1895ca45-3ea2-413b-af87-bc7b06021f5a",
    "execution_millis": 4,
    "execution_start": 1605009995371,
    "output_cleared": false,
    "source_hash": "d112ab9a"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image size =  (111, 126, 126)\n",
      "Image min and max:  0.0 0.94772357\n"
     ]
    }
   ],
   "source": [
    "arr_resampled_dose = sitk.GetArrayFromImage(img_resampled_dose)\n",
    "print('Image size = ', arr_resampled_dose.shape)\n",
    "print('Image min and max: ',  np.amin(arr_resampled_dose), np.amax(arr_resampled_dose))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "cell_id": "00012-f5fc639c-51ac-4afd-a3c1-fcfdb6fd69da",
    "execution_millis": 216,
    "execution_start": 1605009996523,
    "output_cleared": false,
    "source_hash": "2ddd38d8"
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0903a09dbb0443d59d1421c92b848842",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fc6cef4ebbf848c1900360c103ef75ed",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=58, description='nslice', max=111), FloatSlider(value=0.7, description='…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig2, ax2 = plt.subplots()\n",
    "def show_fusion(nslice=58, opacity=0.7):\n",
    "    ax2.imshow(arr_ct[:,:,nslice], cmap=plt.cm.gray)\n",
    "    a = arr_resampled_dose[:,:,nslice]\n",
    "    b = np.ma.masked_where(a <= 0.001, a)\n",
    "    ax2.imshow(b, alpha=opacity, cmap=plt.cm.hot)\n",
    "    \n",
    "interact(show_fusion, nslice=(0,len(arr_ct)), opacity=(0,1,0.1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00013-88abd044-d8f6-4c7d-89ed-5f2f90554e4a"
   },
   "source": [
    "# Part 2: 1D plot from 3D proton dose"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "cell_id": "00014-ffb37061-1c2d-46ab-8948-5c6c9e36f210",
    "execution_millis": 5,
    "execution_start": 1605010005048,
    "output_cleared": false,
    "source_hash": "8ca55087"
   },
   "outputs": [],
   "source": [
    "filename_edep = os.path.join(str(folder), './output_ref/3d-pat-proton-Edep.mhd')\n",
    "img_edep = sitk.ReadImage(filename_edep)\n",
    "arr_edep = sitk.GetArrayFromImage(img_edep)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "cell_id": "00015-f9dd3c8c-a10c-4b90-b6a8-f18f6c5f0ede",
    "execution_millis": 180,
    "execution_start": 1605010005853,
    "output_cleared": false,
    "source_hash": "9b152455"
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ea2d4173a3684c0791f1674c76826776",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Declare a single figure (one row, one column)\n",
    "fig3, ax3 = plt.subplots()\n",
    "moy_dose=np.zeros(arr_dose.shape[0])\n",
    "moy_edep=np.zeros(arr_edep.shape[0])\n",
    "for i in range(arr_dose.shape[0]):\n",
    "    subarr_dose=arr_dose[i,:,:]\n",
    "    moy_dose[i]=np.mean(subarr_dose)\n",
    "    subarr_edep=arr_edep[i,:,:]\n",
    "    moy_edep[i]=np.mean(subarr_edep)\n",
    "\n",
    "# X values from 0 to n\n",
    "# n is the number of slices\n",
    "n = arr_dose.shape[0]\n",
    "x = np.linspace(0, n, n)\n",
    "\n",
    "c1 = ax3.plot(x, moy_dose, 'k.-', label='deposited dose', linewidth=1)\n",
    "ax32 = ax3.twinx()\n",
    "c2 = ax32.plot(x, moy_edep, 'b.-', label='deposited energy', linewidth=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00016-8979e1f6-8d29-4df8-bcc1-3e884650bc7e"
   },
   "source": [
    "# Other examples of exercices\n",
    "\n",
    "- Plot a depth dose profile with interactive line coordinate\n",
    "- Display uncertainty map\n",
    "- Create another simulation with a different beam angle\n",
    "- Sum and display two different dose maps (two different irradiation fields for example)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "deepnote_execution_queue": [],
  "deepnote_notebook_id": "941a92cc-9788-45bd-8e60-5aeb69a51003",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
